Disordered spinor Bose-Hubbard model 
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We study the zero temperature phase diagram of the disordered spin-1 Bose-Hubbard model in a 
2-dimensional square lattice. To this aim, we use a mean field Gutzwiller ansatz and a probabilistic 
mean field perturbation theory. The spin interaction induces two different regimes corresponding 
to a ferromagnetic and antiferromagnetic order, fn the ferromagnetic case, the introduction of 
disorder reproduces analogous features of the disordered scalar Bose-Hubbard model, consisting in 
the formation of a Bose glass phase between Mott insulator lobes. In the antiferromagnetic regime 
the phase diagram differs more from the scalar case. Disorder in the chemical potential can lead to 
the disappearance of Mott insulator lobes with odd integer filling factor and, for sufficiently strong 
spin coupling, to Bose glass of singlets between even filling Mott insulator lobes. Disorder in the 
spinor coupling parameter results in the appearance of a Bose glass phase only between the n and 
n + 1 lobes for n odd. Disorder in the scalar Hubbard interaction inhibits Mott insulator regions 
for occupation larger than a critical value. 

PACS numbers: 64.60.Cn,03.75.Mn,67.85.-d 



I. INTRODUCTION 



Spinor Bose-Hubbard (BH) models describe strongly 
correlated lattice systems where bosons have internal an- 
gular momentum whose orientation in space is not ex- 
ternally constrained. Bosonic interactions are sensitive 
to the spin degree of freedom leading to a rich vari- 
ety of orderings in the ground state at zero tempera- 
ture. In atomic gases, the spin degree of freedom cor- 
responds to the manifold of degenerate -in absence of 
an external magnetic field- Zeeman energy states associ- 
ated to a given hyperfine level F, i.e. {\F,mp)} where 
niF = —F, ..,F. In this context, we identify the spin of 
the atom with the hyperfine quantum number F. Like 
in the scalar case, ultracold atomic spinor interactions 
can be parametrised by two-body short range (s-wave) 
collisions. Due to the rotational symmetry, two-body 
collisions between atoms depend only on their total spin 
and not on its orientation. Moreover, symmetry argu- 
ments impose that the collisions between two identical 
bosons in a hyperfine spin level F are restricted to total 
even spin S = 2F,2F — 2,...,0. Different properties of 
spinor condensates in a single trap has been discussed 
[l|-[j] . The confinement of the particles in a lattice leads 
to an enhancement of the interactions, pushing the sys- 
tem to a strongly correlated regime. As it happens in the 
scalar BH case ^] , the competition between the different 
energy scales present in spinor BH models determines 
the ordering properties -quantum phases- of the ground 
state. Modifying the energy ratio between the hopping 
and interactions allows to cross a quantum phase transi- 
tion between a spinor superfluid (SF) condensate and a 



Mott insulator (MI) state 

The crucial effects of the disorder in condense mat- 
ter systems were advanced in the seminal contribution 
of Anderson [9|, predicting an exponential localization 
of all energy eigenstates of a single particle in a peri- 
odic potential when additional impurities are added to 
it. It took several years to recognize the enormous con- 
sequences Anderson's result had, but nowadays it is well 
established that disorder, and specifically quenched dis- 
order (i.e. frozen during the typical time scales of the 
system), is an essential ingredient in condensed matter 
systems and related topics as conductivity, transport, 
high-Tc superconductivity, neural networks, insulating 
phases or quantum chaos to mention few examples (see 
[lOl and references therein). Disorder is intrinsically dif- 
ficult to treat firstly because, in order to characterize 
the system, one should average over different realizations 
of disorder which is usually a hard task. Secondly, dis- 
ordered systems often develop a complex landscape of 
low energy states making the problem of minimization to 
find ground states very involved. Thirdly, they incorpo- 
rate often fractal and ultrametric structures, all together 
making the problem of simulating quantum disordered 
systems a very complex one. 

In recent years, it has become clear that ultracold 
atoms offer a new paradigm of disordered systems, due 
to the fact that random or quasi random disorder can be 
produced in these systems in a controlled and reproducible 
way. Standard methods to achieve such controlled disor- 
der are the use of speckle patterns [ll|, [l^l which can be 
added to the confining potential, or optical superlattices 
created by the simultaneous presence of optical lattices 
of incommensurate frequencies |13l4l5l |. Other methods 



include using an admixture of different atomic species 
randomly trapped in sites distributed across the sample 
and acting as impurities [la, [l3] , or the use of inhomoge- 
neous magnetic fields which modify randomly, close to a 
Feshbach resonance, the scattering length of the atoms in 
the sample depending on their spatial position [la, U^ ■ 

Strongly correlated bosons in a lattice in the presence 
of external random potentials were first considered in [5] 
where the phase diagram in the i — /i plane of the sys- 
tem, fi being the chemical potential, was worked out. The 
three possible ground states predicted were: (i) an incom- 
pressible MI with a gap for particle-hole excitations; (ii) 
a gapless Bose-glass (BG) insulator with finite compress- 
ibility and exponentially decaying superfluid correlations 
in space; and (iii) a SF phase with the usual off-diagonal 
long range order. Previously, the onset of superfluidity 
in a random potential in ID was studied in [20|, con- 
sidering hard core bosons and using a mean field the- 
ory including quantum fluctuations, and in [2l|, where a 
renormalization group approach was developed to study 
a one-dimensional system of interacting bosons in a ran- 
dom potential. In recent years it has been shown that the 
question of the simultaneous presence of disorder and in- 
teractions constitutes an important and complex many 
body problem that is still far from being well understood 
(for a review see [22|). 

Here we address the effects of disorder in the strongly 
interacting spin-1 BH model in two dimensions (2D). For 
spin-1 systems, the short range two body collisions lead 
to a spin independent effective coupling strength Uq, sim- 
ilar to the scalar case, plus the spinor coupling U2 ■ With 
the help of a Gutzwiller ansatz, supplemented by a per- 
turbative mean field approach, we provide the phase dia- 
gram on different regimes of the phase space determined 
by a spinor coupling and the disorder. The Gutzwiller 
mean field approach is known to give reasonable results 
for the scalar and the spin-1 SF-MI transition in 2D [3]- 
It has also been used to signal in the presence of disorder. 



a BG phase in ultracold scalar bosonic gases [IJ] as well 
as diverse glassy phases in Bose-Fermi mixtures [23]. A 
Gutzwiller mean field approach yield to correct ground 
state for small values of the spinor coupling, since it ne- 
glects correlations between different sites and thus is not 
precise enough in determining accurately the boundaries 
between distinct quantum phases. Nonetheless it pro- 
vides a valuable estimate on the physics of the system 
and permits easily to include the effects of disorder going 
beyond the homogeneous mean field approach. 

One may argue that the mean field approach could 
work even better in the three dimensional case. However, 
the necessarily inhomogeneous, disordered systems are 
then much harder to treat being computationally very 
demanding. For that reason we restrict ourselves to the 
2D case only as in the earlier studies [ij, |23, [2J . 

Our main results can be summarized as follows. In the 
non disordered case, and for U2 > 0, we confirm previous 
findings [s, [3, 123] consisting in: (i) a first (second) order 
phase transition from MI to SF for even (odd) occupa- 



tion numbers in the region U2/U0 < Uc (with Uc — 0.2) 
and t <S^ 1/2, where t denotes the hopping; (ii) a second 
order phase transition from even occupation MI lobes to 
SF if U2/U0 > Uc- In the presence of disorder in the 
chemical potential the above effects, (i) and (ii), persist 
together with the appearance of a BG phase between the 
MI lobes. For U2/U0 > 0.5, odd occupation MI lobes 
disappear while even lobes survive and the correspond- 
ing BG is formed only by singlets between the remain- 
ing lobes. Also disorder can make the odd occupation 
MI lobes to disappear but the BG phase is nematic if 
U2/U0 < 0.5. Assuming disorder in the U2 coupling we 
observe that the BG phase appears only between every 
second pair of lobes and we explain such a peculiar be- 
haviour using perturbation theory in the vanishing tun- 
neling limit. On the other hand, disorder on the spinless 
term of the interaction coupling reproduces qualitatively 
the results found for scalar gases 18]. 

The paper is organized as follows: In section |ll] we 
introduce the spin-1 BH model and shortly review the 
different phases in the homogeneous case (without disor- 
der). In Sec. Ill Al we discuss first the exact phase diagram 
in the absence of tunneling to grasp the features of the MI 
phase. In Sec lHBI we comment the perturbative results 
for small tunneling, while in III CI and IIIDI we derive a 
mean field phase diagram for finite tunneling using both, 
mean field perturbation theory (MFPT) and Gutzwiller 
mean field approach. In section IIIII we analyze in detail 
the effects of disorder. Two types of disorder are consid- 
ered here, disorder on the on-site energies, resulting from 
a random external potential, and disorder on the interac- 
tions both on the scalar and the explicit spin dependent 
part. We calculate the phase diagram in the disordered 
case using both, a Gutzwiller ansatz and MFPT. Finally, 
in Sec. IIVI we present our concluding remarks and open 
questions. 



II. BOSE-HUBBARD MODEL FOR SPIN-1 
BOSONS 



Low energy spin-1 bosons loaded in optical lattices suf- 
ficiently deep so that only the lowest energy band is rel- 
evant can be described by the spinor BH model. The 
corresponding Hamiltonian is [6|: 



{ij),c 
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(1) 



where {i,j) indicates that the sum is restricted to nearest 
neighbors in the lattice and aj^. (Sio-) denotes the creation 
(annihilation) operator of a boson in the lowest Bloch 
band localized on site i with spin component a — 0,±1. 
The first term in ([IJ represents the kinetic energy 
and describes spin symmetric hopping between nearest- 



neighbor sites with site independent tunnehng amph- 
tude t. The second and third term account for spin 
independent and spin dependent on site interactions, 
respectively. These energies at site i are defined as 
Uo^2 = Co, 2 / dr'w'^{r — fi) with cq = ATTh^{ao + 2a2)/3TO 
and C2 = 4Trh^{a2 — ao)/{3m), where as is the s-wave 
scattering length corresponding to the channel with to- 
tal spin 5* [l|, I2I and w{r — fi) is the Wannier function 
of the lowest band at site i. While the second term of 
dD) is spin independent and equivalent to the interaction 
energy for scalar bosons, the third term represents the 
energy associated with spin configurations within lattice 
sites with 



<t<t'=0.±1 



.F^ 
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(2) 



being the spin operator at site i and F the traceless spin-1 
matrices. The explicit form of the spin operator Si reads 



S, 



rii — n_i 



S - — 



fal+O-i 



lAao + H.c. 



—a\ + a_^ I Go — H.c. 



(3) 



S's components obey standard angular momentum com- 
mutation relations [5/, Sj\ = ieijkSk- Note that the spin- 
interaction term favors a configuration with total mag- 
netization zero (denoted as polar and sometimes antifer- 
romagnetic) for U2 > and ferromagnetic for U2 < 
[ij-yl . In the grand canonical approach the total number 
of particles is controlled by the last term of ([T]) where /i 
is the chemical potential and 



E 



(4) 



=o,±i 



is the total number of bosons on site i. Hamiltonian ([T]) 
can be straightforwardly derived from the microscopical 
description of bosonic atoms, with a hyperfine spin F — 
1, loaded in a deep optical lattice and considering the 
two-body short range (s-wave) collisions. More details 
about the derivation can be found in [l|, ^ 0, |2a 123 ■ 

Notice also that, since the orbital part of the wave func- 
tion in one lattice site is the product of Wannier functions 
for all the atoms, it is symmetric under permutation of 
any two atoms. Therefore, the spin part of the wave- 
function should also be symmetric due to Bose statistics. 
This imposes Si + Ui to be even _28] , being Si and rii the 
quantum numbers labelling the eigenvalues of S^ and fii. 

As in the scalar case, the spinor BH system exhibits a 
quantum phase transition between superfluid and insu- 
lating states [a, 0|- In the insulating states, fluctuations 
in the atom number per site are suppressed and virtual 
tunneling gives rise to effective spin exchange interac- 
tions that determine a rich phase diagram in which dif- 
ferent insulating phases differ by their spin correlations. 



The appearance of spin mediated tunneling transitions in 
the optical lattice depends clearly on the ratio between 
the different energy scales appearing on the BH Hamilto- 
nian ([TJ. In alkalins, the scattering lengths are such that 
spin-independent interactions Uq are larger than spin- 
dependent ones U2- In such case the SF-MI transition 
depends mostly on the ratio t/Uo. However, inside the 
insulating regime, the value of U2 plays an important 
role ii U2 > t where it competes with the spin-exchange 
interactions induced by small fluctuations of the particle 
number determining the spin structure. On the contrary, 
if i ^ C/2 tunneling acts similarly for all spin components 
and the gas will behave as a strongly correlated scalar 
gas. 



A. The phase diagram at f = 0. 

To better understand the effects of spin mediated in- 
teractions when disorder is present let us first summarize 
the phase diagram at t — (atomic limit) without disor- 
der. In this limit, the Hamiltonian reduces to the sum of 
independent single-site Hamiltonians Hq = ^^ Ho^i with 



Ho,i = -/i"* + ^ E"'*^^* ^ -"^^ + V E (^« ^ ^'^ 



(5) 



Since 



ni,Sl 



: 0, the eigenstates of the single site 

Hamiltonian can be labeled by three quantum numbers 
\Si,mi;ni), such that: 

H^\Si,mi;ni) = £'o(S'i,nj, ^/q, C/2, ^i) 15*4, m^; n^) (6) 

with: 



EoiSi,n^,Uo,U2,^i) = -p^ni + -Uoni{ni-l) 

+ ]^U2[S^{S, + I) - 2n.,] . (7) 

From Eq.Q, one can deduce the structure of the 
ground state of the insulator phases in the limit i = 0. 
For antiferromagnetic interactions, U2 > 0, the minimum 
energy £^™™ is attained with minimum Si, its specific 
value depending of the number of atoms per site. Thus, 
for even filling factor, the minimum spin is zero and the 
state is described as |0i, 0^; n^) with rii even. This state is 
known as spin singlet insulator |29| . If the atom number 
per site is odd, then the minimum spin per site is one 
and the state reads \li,mi]ni). The chemical potential 
region for which each of the two phases are the ground 
states can be found easily from Eq.©: 
(i) MI with n odd and spin 1 on each lattice site is the 
ground state if £'o(l,n) < _Eo(0, n — 1) and EQ{l,n) < 
^0(0, n + l) i.e. when {n - l)t/o < ^Ji< nllo - 2U2. This 
sets an upper bound on the spin coupling U2/U0 < 0.5 
above which the odd lobes cease to exist [231. 
(ii) MI with n even and spin on each site are ground 
states ii Eo{0,n) < Eo{l,n~l) and Eo{0,n) < Eo{l,n + 



1) leading to (n - 1)Uq - 2U2 < ^J. < uUq for U2/UQ < 
0.5. For higher values, odd lobes do not exist and the 
stability conditions read i?o(0,n) < £'o(0,n — 2) and 
Eo{Q,n) < i?o(0,n + 2). The last two conditions set 
an n-dcpendcnt upper bound on the maximum value of 
U2/U0 < (?i + 1/2). 

The ferromagnetic side of the diagram is easily calcu- 
lated imposing an integer number of particles and realis- 
ing that the minimisation of the energy implies maximum 
spin value i.e. Si = rii. 

The exact phase diagram in the {U2/Uo,fJ./UQ) plane 
is displayed in Fig. [T] providing the width of the MI lobes 
at t = as a function of U2/U0. It is interesting to note 
that the right boundary of even lobes in the range < 
U2/UQ < 0.5 does not change with t/2. This fact leads 
to a stability with respect to disorder in this parameter, 
as we shall show in the Sec. IIII CI corresponding to the 
absence of the BG phase between lobes with occupation 
n and n + 1 with n even in the presence of disorder in 
U2 ■ In the antiferromagnetic region, for f/2 large enough, 
odd lobes disappear while even lobes broaden. In the 
ferromagnetic case the lobes shrink as 1 1/2 1 increases and 
disappear ior U2 = —I- 



We restrict ourselves to 2D systems and revise the re- 
sults from [y]. The MI lobes for odd n are in a nematic 
phase characterized by zero expectation value of all the 
spin components but broken spin symmetry {Sf^') = 
and (Sf^) = (Sfy) = 1 (nemati c p hase for 3D spin 
systems has been studied also in [30, |311)- The corre- 
sponding state is well described by the mean field ansatz 
\ip) — Yli \Si = l,mi = 0). For n even (even MI lobes) 
and sufficiently large U2, the spin-dependent term in 
the Hamiltonian dominates and a singlet configuration 
lip) — J^j l^i — 0,mi = 0) is realized. However, when 
t'^ oc U0U2 tunneling may effectively couple S — and 
S — 2 states leading again to a nematic state. As shown 
in [6\, a first order transition may turn place within the 
MI lobe between the singlet configuration (for low t) and 
the nematic state (at higher t values) with a critical tun- 
nelling rate fulfilling zt'^ = O.5U2U0 for n — 2 where z is 
the number of neighbors. Such phase transition may take 
place only if t^ is sufficiently small so that the MI lobe 
exists at this value, otherwise a singlet MI - SF phase 
transition occurs first and the nematic state may not be 
formed. 




^/Uo 



FIG. 1: Phase diagram of the spinor F — 1 BH model in the 
limit t — 0. Each region corresponds to a MI phase with a 
different occupation number. 



B. Perturbative approach for small t 

For small but finite tunneling t/Uo and \U2\/Uo ^ 1, 
it is possible to perform perturbation theory and derive 
an effective Hamiltonian to second order in t/Ua \o], that 
permits to study insulating phases. The explicit form of 
the effective second order perturbation Hamiltonian de- 
pends on the number of bosons per site, odd or even. A 
mean field theory with a product state ansatz, applied 
to the effective Hamiltonian provides the following char- 
acter of insulating states for the low t limit. For U2 < 
the MI lobe with n bosons is ferromagnetic with S — n. 
The situation is richer for antiferromagnetic ordering and 
depends on the dimension of the system. 



C. Standard Mean Field Perturbative Approach 

The MI-SF transition for spin-1 has also been stud- 
ied using standard mean field perturbative approach by 
Tsuchiya et al. Q. We describe these results in more 
detail since they are a starting point for our study of 
the effects of disorder. Neglecting second order fluctua- 
tions of the bosonic annihilation and creation operators 
we obtain the condition {al ^ — {ai^a)){aj^a — {a-j^a)) — 
which allows to decouple the hopping term as aj^aja — 

i^ia'^jiy + '^LV'icr — i^ia-i^ja, whcrc wc havc introduced the 
superfluid order parameter ipja- — {aja) which, in a ho- 
mogeneous lattice, is site-independent. The Hamiltonian 
reduces now to a sum of local terms Hmf = "Ylii ^i with 



h — —tz 



Uo 






[{i'^al + ip*„aa) - iV'^r 



fj,n 



+ — h{h — 1) + 



El 
2 



2n 



(8) 



where the site index i has been dropped since we are 
considering here an homogeneous system, and z denotes 
the number of nearest neighbours. The superfluid or- 
der parameter tpa has to be determined by minimiz- 
ing the free energy / = — l//31ogTr exp (— /3ft,) , where 

f3 = 1/KbT being Kb the Boltzmann's constant and 
T the temperature. Here, since we are interested only 
at the zero-temperature properties, the former condition 
reduces to the minimization of the ground state energy 
-E'Gs(V'cr) = {GS\h\GS) with the self-consistent condi- 
tion {GS\ flcr \GS) = ■0(T- For sufficiently small t we can 
apply perturbation theory, h = Hq -\- V{t), and use as a 
basis the eigenstates of -ffo ©• The perturbation term 



is given by 



(9) 



Let us focus on the antiferromagnetic case U2 > 0- 
A tedious but straightforward calculation of the matrix 
elements of the perturbation leads to the phase bound- 
aries between the SF phase and the MI phase in the 
{fi/Uo,t/Uo) plane for a given value of the U2/U0 cou- 
pling. Notice that being V ex (a J. + Ha) only even terms 
on the perturbation expansion survive. As derived in [Tj 
the ground-state energy up to second order is for odd 
occupation number given by: 



£;(2)(5=l,„,i,(7o,C/2,M>.) 



= zt 



I — zt 2J aj(n, C/o, C/2,/^) 



Y.\^a\\m 



and for even occupation 



= zt 



zt 



1- Y X! 7j("-, C^o,t^2,M) 



J = 1.2 



J2\i^.\\in) 



with 



ai{n,Uo,U2,fJ-) = 

a2in,Uo,U2,fJ.) = 

a3{n,Uo,U2,n) = 

a4in,Uo,U2,n) = 

ji{n,Uo,U2,fi) = 
j2{n,Uo,U2,fi) = 



n + 2 



4(n- 1) 

l'o5n-l,2-n,l{Uo,U2, HY 

n + 1 

4(n + 4) 

155„+l,2;ri,l(^0,C^2,A«)' 

n + i 

5n+l,l;n,t){Uf), U2, fJ.) 

n 

^n-l,l;n,o{Uo, U2, fJ-) 



(12) 



(13) 



and Sir;nAUa,U2,fJ-) = Eo{l,r,Uo,U2, ti-) - 

Eo{s,n,Uo,U2, n)- Minimisation of the energy for 
a finite order parameter (corresponding to SF) is 
achieved when the expressions on the parenthesis in 
([TU]) and pT|) are negative. On the contrary, the MI 
phase, corresponding to zero order parameter is asso- 
ciated to a positive value of such expressions. Hence, 
the phase boundaries between the SF and the MI in 
the {^/Uo,t/Uo) plane, for a given value of the spin 
interaction U2 are given by: 



1 



odd 



2Ej = 1,4"j("'^0,^2,M) 



(14) 



zJ2i=l,2l3{n,Uo,U2,fi) 



(15) 



Notice also that the dimensionality of the lattice is in- 
cluded through the parameter z which indicated the num- 
ber of nearest neighbours. 

The analysis of the ferromagnetic regime {U2 < 0) can 
be done in the same way imposing the condition S = n. 
Since in this case all the spins are aligned, we can consider 
only one of the components m — ±5 in the perturbative 
expansion. A straightforward calculation leads to the 
following explicit expression for the MI to SF boundary: 

jn + nU2 - ^^) [j-l + n) {1 + U2) - ^J.] 
^'''°~ z(l + t/2 + M) ■ ^ ' 



D. Variational Gutzwiller approach 

The variational Gutzwiller approximation is a non per- 
turbative approach, where the wave function takes a form 
of a product over all M sites of the lattice 

M rimax n S 

lip) =Y[^ gi{n)^fi{S,n) ^ h,{S,m,n)\S,m,n) 



= 1 n=0 



S=0 



(17) 

where gi, hi, fi are the variational coefficients to be de- 
termined by minimizing the BH Hamiltonian ([1} with the 
above ansatz. That implies decoupling in the tunneling 
term (aJ^Ojo-) = n^o-^ -I- (a|^)(a-,CT)(l - S^j). Observe 
that for consistency of notation we should rather use ipia- 
instead of {aja). The Gutzwiller variational state is a 
product state of on-site wave functions so it cannot re- 
produce intersite correlations or entanglement between 
different sites. Being a generalisation of the standard 
mean field approximation, the Gutzwiller ansatz is ex- 
pected to be exact in the limit of infinite dimensions. To 
mark the limits between the SF and MI phases in the 
Gutzwiller approach, we recall that the MI phase pre- 
vails for small hopping amplitude and it is characterized 
by a finite gap in the spectrum and zero compressibility 
defined as k = q^ with 



(18) 



and N the total number of bosons. On the contrary, 
in the SF phase bosons are delocalized and a current 
flow is possible. This phase is characterized by a flnite 
compressibility, gapless excitations and off-diagonal long 
range order accompanied by a non vanishing order pa- 
rameter. Since the order parameter is not directly mea- 
surable, it is important to deflne experimental observable 
quantities marking the SF phase. These are typically the 
superfiuid fraction ps and the condensate fraction pc 
[13| (although it has been proposed recently the mea- 
sure of the compressibility directly [S^l). The superfiuid 




fraction can be evaluated imposing a phase gradient in 
the tunnehng corresponding to a current flow while the 
condensate fraction is defined as the hi ghe st eigenvalue 
of the one particle density matrix [ij, |lj]. In an ho- 
mogeneous case the site dependence can be omitted. It 
is measured experimentally by means of an interference 
density pattern giving coherent peaks in the SF phase 
[33| . Notice that in the mean field approach, and so in 
the Gutzwiller ansatz, the condensed fraction decouples 
and both quantities, superfiuid fraction and condensed 
fraction are related to the average ipa — (o-aj), that can 
be taken as the order parameter, its value being zero in 
the MI phase and finite in the SF. 

In Figure [5] we display the pc calculated both with 
the Gutzwiller ansatz and with the perturbation mean 
field approach (|14m5p (solid line) for different values of 
the parameter U2/U0 (left column). Observe the different 
behavior between odd and even lobes (as described in the 
previous section). With increasing values of U2, the even 
lobes start to dominate while the odd lobes shrink. We 
have checked numerically that for U2 = 0.5Uo the odd 
MI lobes disappear, in agreement with the t — pre- 
dictions of the previous section. For small U2/U0 ratios, 
there exist a discrepancy between the perturbative mean 
field and Gutzwiller predictions for the boundaries of the 
even lobes already reported in [3J|. This discrepancy is 
correlated with the character of MI-SF transition as visu- 
alized in the right column of Fig. [51 where the condensate 
fraction is shown for selected /i = const lines correspond- 
ing to the tips of the lobes in the corresponding phase 
diagram. 

For U2/U0 < 0.1 (Fig.[2]first and second row) the con- 
densate fraction is continuous across the phase transi- 
tion for odd lobes (corresponding to second order phase 
transition) while it reveals a discontinuous jump, charac- 
teristic of the first order phase transition for even lobes. 
For U2/U0 > 0.3 (Fig. [2] bottom row) only second or- 
der SF-MI transitions from both odd and even lobes are 
observed. In between these values it is not easy to char- 
acterize the order of the phase transition. An exhaustive 
numerical analysis shows that the value where the tran- 
sition passes from first to second order is approximately 
U2/U0 = Mc - 0.2. 

The observation of a first order phase transition in the 
even lobes - where MI is formed by singlets on each site- 
is not new and has been also pointed out in the mean field 
analysis of [3, l3J| in 2D, as well as in Quantum Monte 
Carlo (QMC) calculations ^35] in ID. Recall also that 
another first order transition between singlet and nematic 
phases has been predicted within the MI lobes 6] by 
using a restricted MF ansatz in the effective perturbative 
Hamiltonian. 

Looking at the state provided by the Gutzwiller ansatz 
near the even lobes' tips (Fig. [3]), we observe that in the 
MI only the S = component is relevant, while in the SF 
phase the state becomes a linear combination oi f(S = 
2,n)\S = 2,m, = 0, n)+f{S = 0, n) jS" = 0, m = 0, n) with 
|/(S' = 2,n)|2 + |/(S' = 0,n)|2 close to 1. A close inspec- 



tion of the coefficients of the Gutzwiller ansatz (ITT)) shows 
that, for U2/U0 < 0.1, f{S — 2,n) assumes a finite value 
abruptly (see Fig|3](a) (b) (c)). Both states tend to con- 
tribute equally in the limit case in which U2 = (scalar 
case). Notice that this is also the origin of the discrep- 
ancy with the MFPT result where only states with S ^ 
are taken into account in the energy corrections. Even if 
we find that, in the MI, the state is singlet, the second or- 
der phase transition indicates a metastability inside the 
lobe of a nematic phase. A rough explication of this effect 
can be made for the lobe corresponding to n — 2, noticing 
that a configuration with S = 2 starts to become favor- 
able when the kinetic energy becomes comparable with 
Eo{S = 2,n) - Eo{S = 0,n) = 3U2, so for zti ~ 6U2. 
This value can be compared with the MI tips obtained in 
the MFPT zt2 = (Uq + 2U2) [{2n -h 3) - V4n2 -|- 12nJ . If 
ti < t2 the kinetic energy reduces the lobe with respect 
to the MFPT prediction and, as soon as the metastable 
nematic state becomes stable, SF phase appears discon- 
tinuously. On the other hand, if ^2 < ^i the system be- 
comes SF before the appearance of the S' = 2 contribu- 
tion and the MI-SF transition is smooth. It is easy to 
verify that ii ~ ^2 for U2/U0 ~ 0.145 which is not too 
far from the Uc ~ 0.2 mentioned above. In contrast, for 
U2/U0 > Uc Gutzwiller and MFPT approaches coincide 
and effectively the contribution of the state S' = 2 is irrel- 
evant close to the tip of the lobe (see Figl3](d)) and the 
transition rather than a phase crossing becomes again a 
second order phase transition. 



III. DISORDER IN SPINOR BOSE-HUBBARD 
MODEL 



As discussed in the introduction, the presence of disor- 
der in the BH model allows, apart from MI, for another 
insulating phase, the BG phase. The characteristics of 
the phase diagram depend on the way the disorder is in- 
troduced. Here we will study the effect of two different 
kinds of disorder: disorder in fj, (Sec. IIIIB[) . and disorder 
in the interactions U2 and Uq (Sec. IIIIC[) . 

Diagonal disorder can be taken into account by adding 
to the Hamiltonian ([IJ a local term such 



Hd — H + 2_^ Hdis (ci) 



(19) 



where e^ is a random variable defined for every site i with 
a given probability distribution p(e). While, depending 
on the origin of the disorder, different p{e) may be con- 
sidered (see e.g. [321 )■ Here we consider the simplest 
uniform distribution with — A < e^ < A and the cases in 
which the disorder is equivalent to add a random term 
to one of the variables ^, U2 or Uq. 

Notice that the addition of a site dependent disorder 
introduces inhomogeneity into the system. Thus, neither 
the mean field nor the Gutzwiller ansatz reduce to a "sin- 
gle site" effective Hamiltonian. Instead, the mean fields 
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FIG. 2: (Color online). Left panels depict the condensate fraction pc obtained numerically by the Gutzwiller ansatz for, 
(a)C/2/f/o = 0.02, {h)U2/Uo = 0.1 and {c)U2/Uo = 0.3, in the homogeneous case without disorder, where MI lobes correspond 
to vanishing pc (orange areas). The lobes are compared with the boundaries obtained with the MFPT (solid lines). In the 
right panels is depicted pc as a function of zt/Uo for values of n/Uo corresponding to the lobes' tips. In the transition between 
the MI and SF on the tip one can observe a first order transition for the even occupation lobes in panels (d) and (e) (abrupt 
jump on the condensate fraction) while for lobes corresponding to odd occupation the transition is always of the second order. 
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FIG. 3: Coefficients \f{S = 0, n = 2)f (triangles) and |/(S = 
2, n = 2)1 (squares) of the Gutzwiller state (|17[l as a function 
of zt/Uo. The value p, corresponds to the n = 2 MI lobe's 
tip. The panels refer to different values of spin interaction: 
U2/U0 = 0.01(a), 0.02(b), 0.1(c) and 0.3(d). 



as well as Gutzwiller wave function coefficients become 
explicitly site dependent. 

A Stochastic Mean Field Theory (SMFT), taking into 



account the inhomogeneity of ipi^, has been proposed in 
(36j for the scalar BH. Here we present a more simple 
MP theory, being a limiting case of SMFT, and compare 
it with the phase diagram obtained with the Gutzwiller 
ansatz. Good agreement for the MI boundaries has been 
found, as in the non disordered case, while the BG can 
be seen only by the Gutzwiller ansatz. 



A. Probabilistic Mean Field approach 

A first estimation of the MI lobes in the presence of 
disorder can be obtained by a MFPT, as described in 
the previous section. The generalization to the disor- 
dered case is not straightforward since, as we mentioned, 
the translational invariance is broken and the order pa- 
rameter should be associated to a random variable V'jo- 
defined for each site, with a certain probability distribu- 



tion P{ipj 



Nevertheless, since the disorder is assumed 



to be homogeneous on the lattice, we can introduce a 
simplified MFPT theory taking an average order param- 
eter tpa — J dipjcrPii^ja-)ipjtT- In doing so, we are ne- 
glecting the classical fluctuations of the order parameter 
induced by the d isorder. The self-consistent condition 
reads ij^a- — (o^.a); where the overbar indicates an ensem- 
ble average over the lattice. It is also equivalent, due 
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to the self- aver aging properties of the systera, to an av- 
erage over the random distribution. So the structure of 
the single-site mean field Hamiltonian remains the same, 
providing that one of the parameters changes according 
to V ^ V + ej, where, for a diagonal disorder, v can be 
/i, U2 or Uo- 

The minimization of the average ground-state energy 
which, up to second order corrections, reads 



E{s,n,t,UQ,U2,tJ-,'>jJcr) 



Eo{s,n,Uo,U2,^i) (20) 



determines if ipa- is finite or zero. Notice that, since we 
are neglecting the fluctuations on the order parameter, 
a vanishing ip^ always corresponds to a MI phase, while 
BG cannot be detected, since it has -ipa = but finite 
fluctuations. A more complete analysis needs a more 
complex theory, such as the SMFT [3a, |33| where fluctu- 
ations are taken into account and P{'4>ja) is determined 
self-consistently. 

In our simplified stochastic approach the MI-SF 
boundary is calculated using the equations for the MI 
boundaries (fTH) and (|15p but using the averaged values 
cij and jj instead of aj and jj . 



B. Disorder in /i 

We consider first the disorder in the chemical poten- 
tial corresponding to Hdis{^i) = uf^-i- To study the phase 
diagram, we use the Gutzwiller approach with a lattice 
large enough that self-averaging over the possible dis- 
order realizations is already realized. Now we have to 
distinguish between three phases: SF, MI and BG. As 
before the (disorder averaged) condensate fraction helps 
to find the border between SF and insulator (BG, MI) 
phases. 

For MI, as mentioned in the previous section, both 
the compressibility and the fluctuations in the average 
occupation number vanish within the Gutzwiller ansatz 
approach. The latter simply because the MI is realized 
as a Fock state with the same occupation at each site. In 
a BG phase the wavefunction is again a product of Fock 
states at each site but with different occupations (due 
to local action of the disorder). Thus for a BG, fluctua- 
tions in the average (over sites) occupation number are 
significant. We have checked that we obtain numerically 
practically the same border between MI and BG using 
fluctuations in the average occupation number or by di- 
rectly calculating the compressibility from its definition 
(see the previous section). 

Let us mention that situation is so simple and unam- 
biguous in Gutzwiller approximation only. For finite tun- 
neling the real MI state is not a Fock state and fluctu- 
ations of on site occupation change smoothly across the 
MI-SF transition (see e.g. [33 )■ In experimental situ- 
ation, in addition, atoms are held in an additional trap 
so the density of atoms depends on the position in the 



trap. Then, however, one can use directly compressibility 
measurements for finding MI borders as experimen tally 
shown for fermions [33] and also proposed for bosons [S^] ■ 
Standard time of flight interference patterns then allow 
to determine the condensate fraction. 

Fig. m shows the results obtained for a fixed ampli- 
tude of the disorder, and different values of C/2. Mott 
Insulator lobes correspond to vanishing density fluctua- 
tions and zero compressibility as shown in the left panels 
(a-c). The results obtained from MFPT are also dis- 
played in the panels as a solid line for comparison. As 
in the scalar bosonic case, disorder slightly shrinks and 
separates the Mott lobes and a BG phase appears be- 
tween them. The regions in the (/z/C/g, t/C^o) plane as- 
sociated with the BG phase are obtained by contrasting 
results obtained from the condensate fraction {pc) with 
the zero-density fluctuation regions (MI lobes). The re- 
gions associated to bose glass phase correspond to those 
regions where fluctuations are different from zero (com- 
pressible) but have vanishing condensate fraction. These 
regions are depicted in Fig. 2) (panels d-f ) . In these re- 
gions the single site superfluid parameter can be different 
from zero but has to vanish on average so to destroy the 
off-diagonal interference terms of the pc- 

As one can see in Fig. |4| , no BG appears close to 
the tip of a given lobe, yielding a direct SF-MI transition 
even in the presence of disorder. This is a limitation of 
the mean field approach. Recently it has been claimed 
by means of the noninclusion theorem and supported by 
QMC calculations, that BG always separates SF from MI 
phase [401 . 

The MI phase, in the scalar BH model, disappears com- 
pletely for A > 0.5C/o- This may be easily understood 
from the fact that the maximal possible gap separating 
the ground state and first excited states is, in a homo- 
geneous case and in i — >■ limit, equal to Uq. Thus 
disorder spanning [— C/o/2, C/o/2] interval effectively fills 
up the gap, producing a disordered gapless medium 41^. 
The same argument may be used for odd and even lobes 
in the spinor case. For the even lobes the maximal gap 
is C/q -I- 2C/2 when U2 < 0.5Uo while for odd occupation 
lobes the maximal gap is Uq — 2U2- Thus the critical dis- 
order for the disappearance of the odd occupation lobes 
is Ao = Uo/2 — 1/2- In Fig. Hjthe values of critical dis- 
order are, from top to bottom, Aq — 0.48?7o, OAUq and 
0.2C/o- Then, A > Aq only for the last case, where we 
see the disappearance of the odd occupation lobes. It is 
interesting to note that when the odd filling MI is sup- 
pressed, the BG is nematic for U2/U0 < 0.5 while it is 
formed by singlets for U2/U0 > 0.5. This fact can be seen 

in Fig. [S] where the averaged ( S^ ) is plotted as a func- 
tion of /i/C/oi for zt/Uo = 0.02 and four different values 
of U2. Comparing this plot with FigHl one can see that, 
for U2/U0 < 0.5, MI phases correspond to constant value 

of S, either (sA = 2 (odd lobes) or /s^N = (even 

lobes) . Outside this constant values the associated phase 
is BG. For U2/U0 = 0.3, where odd filling lobes exist in 
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FIG. 4: (Color online). Left column panels report the average density fluctuations v n^ 
approach. MI lobes, corresponding to vanishing fluctuations (orange areas), are compared with the probabilistic mean field 
prediction (solid lines). Right column panels show the corresponding condensate fi-action in comparison with the Gutzwiller 
MI lobes (solid lines). The zero-condensate fraction areas (orange areas) outside the MI lobes correspond to BG phase. For all 
panels, random disorder in the chemical potential with A = O.SUo is considered. The different panels correspond U2/U0 = 0.02 
(a-d), U2/U0 = 0.1 (b-e) and U2/U0 ~ 0.3 (c-f). Observe the disappearance of the odd filling MI lobes for the largest U2/U0 
ratio in agreement with the simple estimate given in text. 



the ordered case but are suppressed by the disorder, the 
BG has 0<S'<1(0<(S2)<2) corresponding to a 



nematic phase (since (S'^) =0 and (S^^ ^ 0). On the 

other hand, for U2/U0 > 0.5 both, the even filling MI and 

the BG, have (S'^) = 0, meaning that they are formed 

by singlets. So disorder can destroy insulator with odd 
fiUing, but only for U2/U0 > 0.5 singlet Bose Glass is 
formed. 

Notice that, as in the non-disordered case, the 
Gutzwiller ansatz closely coincides with the SMFA for 
the boundaries of the odd occupation lobes while it dis- 
agrees for the even ones for sufhciently small U2 (FigH]). 
Based on the intuition obtained from the case without 
disorder we may again associate the disagreement with 
the hidden first order transition. Such a situation occurs 
for t < y/lhXh). For larger U2 MFPT and the Gutzwiller 
approach produce practically identical results. 
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FIG. 5: Total spin average (o) as a function of h/Uq with 
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C. Disorder in U2 and Uo 

One could imagine that disorder in the on-site interac- 
tions U2 and Uq can be experimentally realized, in prin- 



ciple, using optical Feshbach resonances [19|, |42h44| | (the 
application of magnetic field in a standard Feshbach res- 
onance technique would additionally modify the system 
due to e.g. Zeeman level splitting). However the opti- 
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cal Feshbach resonance introduces losses dues to sponta- 
neous emission from the intermediate state [43,|4l] so it is 
not clear at all whether the timescale for losses would al- 
low for realizing the ground state of the system. Very re- 
cently, however, another microwave-Feshbach resonance 
technique has been suggested [431 ■ This method uses res- 
onant microwave driving between ground state sublevels 
to tune the scattering length. Since excited states are 
not involved in this method no additional losses due to 
spontaneous emission are expected. Although this ap- 
proach is up to now a theoretical proposal, it seems to 
be a promising candidate for tuning the interactions in a 
stable way without the application of the magnetic field. 

A small local fluctuation in the laser tuning (assum- 
ing optical scheme with the reservation discussed above) 
or the microwave tuning (in the method of [45|) Suj in- 
troduces fluctuations in Uq and U2 so, in principle, dis- 
order should be considered in both parameters. Since 
the variations 6as of the two scattering lengths are cor- 
related variables (both being function of Suj), we could 
manage to compensate them so to have an almost vanish- 
ing sum or difference. If, for instance, we put the system 
between the two Feshbach resonances, a small detuning 
will increase one scattering length and decrease the other 
one. So, if the condition Sao{6uj) + 26a2{Suj) ~ holds, 
only disorder in U2 can be considered, on the contrary, if 
6a2{Suj) — SaQ{Suj) ~ we can consider only disorder in 
Uo- 

Let us start considering disorder in 1/2, so assuming 
that for each site 1/2 = 1/2 + e^ where Ci takes a random 
value in the interval [—A, A]. Throughout this section 
we assume that A < |f/2| so to consider all U2 of the 
same sign, negative or positive, for the ferromagnetic or 
the antiferromagnetic cases, respectively. 

Figure ini shows the effect of disorder for U2/U0 = ±0.1 
and A/Uo = 0.06. In the ferromagnetic case (plots (a) 
and (c)), disorder in U2 has the same effect as the disor- 
der in ^. Similar to the scalar case, MI lobes are shrunk 
and BG phases appear between them. In contrast new 
features emerge in the antiferromagnetic case (plots (b) 
and (d)), where BG is formed only between lobes corre- 
sponding to n and n+l occupations with n-odd. No effect 
of disorder is visible between n and n + l MI lobes for 
n-even. A very simple explanation of that behavior may 
be obtained from Fig.[Tl Note that for U2 € [0, C/o/2] the 
border separating n and n + l occupations for n-even in 
the H — U2 plot is vertical (in i = limit). Thus changes 
(e.g. fluctuations) in U2 do not modify the chemical po- 
tential at which the density changes. For n odd in this 
range of U2 the border is tilted, thus fluctuations in U2 
for a flxed /i change the density value which is favored 
for the ground state. Then, depending on the particu- 
lar value of U2 at a given site the density for the ground 
state changes. Interestingly, this picture, established for 
i = 0, seems to hold also for finite t as no BG is observed 
between the odd and even lobes. Further inspection of 
Fig. [T] reveals that in other possible ranges of U2 the lines 
separating different densities are always tilted - indicat- 



ing possibility of BG creation between the MI lobes. In- 
cidentally, we can also interpret the same figure assuming 
fixed U2 and fluctuating ^ as the case discussed earlier 
in this paper. There are no horizontal lines in Fig. [T] 
thus all density borders are vulnerable to fluctuations in 
/i. This is again consistent with the observation that for 
disorder in /i BG appears between all lobes. 

Finally, in Fig. [7] we show the result obtained for the 
disorder in Uq. As in the previous case we take Uq = Uq + 
€■1 with €i G [—A, A]. The plots report the case with zero 
and finite value of U2 and A/Uq = 0.25. As explained 
in [18], in the U2 = case, lobes with occupation n > 
{1+A/Uq)/(2A/Uo) disappear while the first one remains 
always stable. In our analysis we recover this behavior 
even for finite 1/2- In both cases lobes separate and BG 
appears in between. 



IV. SUMMARY-OPEN QUESTIONS 

We have analyzed the effects of disorder in the spin- 
1 BH model in which the spin interaction induces two 
different regimes, corresponding to a ferromagnetic and 
antiferromagnetic order, focusing mainly on the antifer- 
romagnetic case, where the phase diagram differs more 
from the scalar case. We have considered here both, dis- 
order introduced to the chemical potential (correspond- 
ing to an offset of energies at different sites) as well as 
disorder in the atom-atom interactions. As for the scalar 
bosons, we have observed the appearance of a compress- 
ible insulator - the BG phase - its character depending on 
the U2/U0 ratio. For small U2 when Mott states with an 
odd number of atoms per site (also termed nematic since 
they have the mean value of all components of the spin 
equal zero, but a non vanishing singlet projection)exist 
in the absence of disorder, we expect the BG to be also 
nematic. For large U2 however, when odd MI lobes do 
not exist already in the absence of disorder, we find a BG 
of singlets, a novel phase peculiar to bosons with spin. 

Interestingly enough, in the presence of disorder on 
spinor coupling U2, the system shows robustness against 
BG creation which does not emerge between n and n+l 
MI lobes for n- even. This is traced back to the insensitiv- 
ity of the MI borders to changes in U2 in the /x/C/o— C/2/C/0 
plane observed in the vanishing tunneling limit. 

This work is only the first step towards understanding 
disorder on lattice spinorial bosons. For ID systems den- 
sity matrix renormalization group (DMRG) or its vari- 
ants may be used to go beyond the mean field; work in 
this direction is in progress. For 2D, similar studies may 
be undertaken within QMC. 

Finally, we remark the suitability of these systems 
for spin-glass studies. Notice that if one induces disor- 
der in the U2 coupling not preserving the ferromagnetic 
and anti-ferromagnetic character of the two-body inter- 
actions, i.e. if A > IC/2I, a situation resembling frus- 
tration will appear in this model with antiferro and ferro 
sites randomly distributed along the lattice. Last but not 
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FIG. 6: (Color online). Density fluctuations (left panels) and pc (right panels) for IJi — ±0.1[/o and disorder in U2 A/Uo = 0.06. 
MI lobes compared with the MF results (solid lines). Vanishing pc outside the MI lobes (solid lines), corresponds to the BG 
phase. Panels (a) and (c) correspond to the ferromagnetic case f/2 = — O.lf/o (a-c). The case C/2 = O.lUo is reported in panels 
(b) and (d). 
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FIG. 7: (Color online). Density fluctuations (left panels) and pc (right panels) for disorder in Uo with A/Uo = 0.25 . MI lobes 
compared with the MF results (solid lines). Vanishing pc outside the MI lobes (solid lines), corresponds to the BG phase. 
Panels (a) and (c) correspond to U2/U0 — 0.0, (b) and (d) to U2/U0 = 0.1. 
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least, a more realistic model of fluctuations in the inter- 
actions should be considered taking the details of optical 
Feshbach resonance into account. 
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